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ABSTRACT 



We predict theoretically the nondiffractive propagation of sonic waves in periodic acoustic 
media (sonic crystals), by expansion into a set of plane waves (Bloch mode expansion), 
and by finite difference time domain calculations of finite beams. We also give analytical 
evaluations of the parameters for nondiffractive propagation, as well as the minimum size of 
the nondiffractively propagating acoustic beams. 



PACS numbers: 43.20.Fn, 43.35.Cg, 43.20.Wd 
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I. INTRODUCTION 



The study of the dynamics of waves, even in simple hnear media, initiated hundreds of 
year ago, ever and ever leads to surprisingly new results and insights. One of such "surprises" 
was the discovery of band gaps in the propagation of light in materials with the refraction 
index periodically modulated on the scale of the optical wavelength, the so called photonic 
crystals^. The theory of wave propagation in periodic materials was developed long time ago 
by Bloch and Floquet, and it found many applications in solid state physics, in particular in 
the studies of electronic properties of semiconductors (calculation of valence and conduction 
bands, etc). Nevertheless, the advent of the photonic crystals initiated a revival of the theory 
of wave propagation in periodic media. The creation and control of photonic band gaps^, 
the slowing down of light^, and the photonic crystal waveguides are the main applications to 
the date. Most of these studies concern the propagation of plane waves (not the beams), and 
results in the modification of the temporal dispersion relation (frequency versus propagation 
wavenumber). Later, the strong analogies between the propagation of light and sound (which 
obey similar wave equations) motivated the study of sound propagation in periodic acoustic 
media, the so called sonic or phononic crystals (SC). Many of the results obtained in the 
photonic case have been reported in the sonic case. For a review on this topic, see e.g. Ref. 
4. 

Most of the studies reported above concern the one-dimensional (ID) periodic struc- 
tures, as the ID case, being relatively simple, allows an analytical treatment. The multi- 
dimensional cases (the 2D case as in our present study, or even the 3D case) are much more 
difficult to be accessed analytically. The majority of these studies in multi-dimensional 
case are numeric, as using plane-wave expansion, or finite difference time domain (FDTD) 
schemes. These studies also mostly concern the modification of the temporal dispersion 
characteristics. 

It comes out recently, that the spatial periodicity can affect not only temporal dispersion, 
but also the spatial one, i.e. the dependence of the longitudinal component of the propaga- 
tion constant versus the transverse component. These results (again predominantly numeric) 
lead to so called management of spatial dispersion, i.e. to the management of diffraction 

properties of narrow beams. This idea led to the prediction of the negative diffraction of 
light beams in photonic crystals^, of sound beams in sonic crystals®, and of coherent atomic 
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ensembles in Bose- Einstein condensates in periodic potentials^. In particular it has been 
found recently that between the normal diffraction and negative diffraction regimes a strong 
reduction of the diffraction can be achieved, leading to the so called self-coUimating, or 
nondifTractive light beams^. 

The geometrical interpretation of wave diffraction is as follows: wave beams of arbitrary 
shape can be Fourier decomposed into plane waves, which in propagation acquire phase 
shifts depending on their propagation angles. This dephasing of the plane wave components 
results in a diffractive broadening of the beams. Fig. 1(a) illustrates normal diffraction in 
propagation through an homogeneous material, where the longitudinal component of the 



wavevector depends trivially on the propagation angle, k\\ = /c^ = y |k| — |k^| , where 
k_|_ = {kx,ky). In general, the normal or positive diffraction means that the surfaces of 
constant frequency are concave in the wavevector domain k — {kx, ky, kz), as illustrated in 
Fig. 1(a). The negative diffraction, as illustrated in Fig. 1(b), geometrically means that the 
surfaces of constant frequency are convex in wavevector domain. The intermediate case of 
the vanishing diffraction is illustrated in Fig. 1(c), where the zero diffraction is supposed to 
occur at a particular point in the wavevector domain where the curvature of the surfaces of 
constant frequency becomes exactly zero. Zero diffraction physically means that beams of 
arbitrary width can propagate without diffractive broadening or, equivalently, that arbitrary 
wave structures can propagate without diffractive "smearing" . 

The present study concerns the nondiffractive propagation of sound in periodic acoustic 
materials (sonic crystals). We found, by applying the plane- wave expansion method, the 
existence of nondiffractive regimes similar to those in optics, or to those expected from 
Fig. 1(c). We check the nondiffractive propagation by integrating the wave equations by 
means of the FDTD technique. Moreover, we also present the analytical treatment of the 
problem, leading to analytic relations, which among other are useful for the planning of the 
corresponding experiment, and for designing the possible applications. 

In Section II of the article the propagation of sound is analysed by plane wave expan- 
sion, leading to the spatial dispersion curves, and in particular resulting into the straight 
(nondiffractive) segments of the spatial dispersion curves. In this way the nondiffractive 
propagation regimes are predicted. In the next Section III the FDTD calculations are per- 
formed in the predicted nondiffractive regimes, and the nonspreading propagation of narrow 
beams is demonstrated. Section IV is devoted to the analytical study, to the derivation 
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of analytical relations between parameters for the nondiffr active propagation. Last Section 
contains the concluding remarks, where the results are summarized and also the minimal 
size of the beams propagating nondiffractively is evaluated. 



II. DISPERSION IN SONIC CRYSTALS 

The propagation of sonic waves is determined by the following linear system equations 
P-Q^ = -Vp, (la) 

where B{r) is the bulk modulus, p(r) is the density (both dependent in space), p{r,t) is the 
pressure (which are scalar fields), and v(r, t) is the velocity vector field. 

We define the relative values of the bulk modulus B{r) = B{j)/ Bh and the density p(r) = 
p(r) / p/j, normalizing to the corresponding parameters in the host medium. Then, eliminating 
the velocity field in Eqs. (0), we obtain a wave equation describing the propagation of sound 
in the inhomogeneous medium, 

ii;y%^--(^-^<-')-- 

where r = Cht is a normalized time, that makes the velocity of sound in the host medium 



Ch equal to unity, where Ch = yBh/ph- 

We consider sound beams with harmonic temporal dependence. Then, the steadily os- 
cillating solution has the form p{r,t) = p{r)e^'^'^, which substituted in Eq. ^ leads to the 
eigenvalue equation 

^pW + v(-±^Vp(r))=0^ (3) 

For the subsequent analysis we consider a concrete geometry, where acoustic waves prop- 
agate in a two-dimensional medium, formed by an squared array of solid cylinders, with 
axis along y direction and radius ro, in a host fluid medium. The coordinate r in Eq. Q 
depends now on longitudinal (z) and transverse (x) directions, and V = {d / dx,d / dz) . 

The lattice defined though the centers of cylinders is given by the set R = 
{R = riiSLi + n2SL2; ni,n2 G A^} of two-dimensional lattice vectors R that are generated by 
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the primitive translations ai and a2. The corresponding reciprocal lattice is defined though 
G = {G : G ■ R = 21171] neN}. 

A possible way of solving Eq. is by means of the plane wave expansion (PWE) method, 
which converts the differential equation into an infinite matrix eigenvalue problem, which is 
then truncated and solved numerically. By solving this eigenvalue problem the frequencies 
corresponding to each Bloch wave can be obtained, providing the dispersion relationship 
and band structure of the periodic medium. 

The bulk modulus and density are periodic functions with the periodicity of the lattice, 
and therefore contain all the information of the phononic crystal. This implies that the 
material parameters can be represented by their Fourier expansions on the basis of the 
reciprocal lattice, 

PW-^ = $^PG'e^^-, (4) 

G 



m-' = Y.^G^'''''- (5) 

G 

On the other hand, the solutions p(r) of Eq. (jS)) must be also periodic with the periodicity 
of the lattice (Bloch-Floquet theorem), and can be expanded as 

p(r) = e^''-^5^Pk,Ge^''^ (6) 

G 

where k is a two-dimensional Bloch vector restricted to the first Brillouin zone, and G 
denotes a reciprocal lattice vector. For a square lattice, G = (27r/a)(niei + ^262) with rii 
and n2 integers and a being the lattice constant (minimal distance between the centers of 
neighbor cylinders). 

The coefficients in expansions Q and (jSj) can be obtained from the inverse Fourier trans- 
form. For the inverse of mass density the coefficients result^ 

Po=-JI ^,dv=P^f+{l-f), for G = 0, (7) 
J J P(r) pc 

uc 

which represents the average value of the density, and 

1 r r „iG r / „ 

-1 ^11^ J„ I P^^ 



UC 



Jl 


(|G| 






|G| 





where the integration extends over the two-dimensional unit cell, Ji{x) is the Bessel function 
of the first kind and / = 7r(ro/a)^ is the filling fraction. Exactly the same expressions, follow 
for the coefficients of bulk modulus 6q^, since the expansion has an analogous form. 
In terms of the coefficients of the previous expansions, Eq. (jH)) becomes 

yb^'-o' - Pg-c (k + G) • (k + G')] PC = 0. (9) 

G' 

Equation Q has been numerically solved considering 361 plane waves in the expansion. 
The number of plane waves has been chosen in order to ensure the convergence. Figure 

2 shows the band structure for a square lattice of steel cylinders (pc = 7.8 10^ Kg m~^. 
Be = 160 10^ N m-2) immersed in water {ph = 10^ Kg m'^ Bh = 2.2 10^ N m'^). The 
dimensionless (reduced) frequency Q = u!a/27cch is plotted in terms of the dimensionless 
wavenumber of Bloch vector K = ka/27r. 

From the solutions of Eq. Q we can also compute the isofrequency contours. In Fig. 

3 the results for the first and second bands are shown. In both cases, the curves show 
a transition from convex to concave at a particular frequency. The isofrequency contours 
at the transition point acquire, as shown in the figure, the form of squares with rounded 
corners. Consequently, there exist locally fiat segments of the curve, where, in other words, 
the spreading of the beam will be counteracted by the crystal anisotropy. Similarly as 
for photonic crystals in optics the nondiffractive propagation occurs along the diagonals of 
squares in the first propagation band, and along the principal vectors of the lattices in the 
second band. The "rounded nondiffractive square" is situated around the corner of Brillouin 
zone (denoted by M) for the first band, and around the centre of Brillouin zone (denoted by 
F) in the second band. 

III. NUMERICAL RESULTS 

In order to prove the nondiffractive propagation of sound in the sonic crystal, Eqs. ([Q) 
have been numerically integrated using the Finite Difference in Time Domain (FDTD) 
method. FDTD is nowadays a standard technique^° for solving partial differential equations 
by integrating in time, and by approximating the spatial derivatives by finite differences. 
The incident acoustic beam has been simulated by a plane vibrating surface radiating a har- 
monic wave with variable frequency u, describing an acoustic transducer with a diameter 
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of 3 cm. The front face of the crystal is located close to the source, where the wavefront is 
nearly plane. The medium parameters were chosen to correspond to numerically evaluated 
zero diffraction point (by inspecting the isofrequency curves) in the previous section. For 
the first band [Fig. 3(a)] the isofrequency curve becomes locally fiat for Q ^ 0.54, which 
corresponds to a real frequency of / = Qch/a ~ 154 KHz, and for an incidence along [1,1] 
direction. Under these conditions, the nondiffr active propagation is predicted to occur. The 
result of the numerical simulation for these parameters is shown in Fig. 4 (left column). As 
expected, the beam propagates through the crystal without a visible divergence. For the 
second band, the theory predicts that the frequency for nondiffractive propagation increases 
roughly by the factor of \/2 with respect to the first band, and then occurs at / ^ 217 KHz. 

We note that whereas the beam in homogeneous media broadens sensibly over the propa- 
gation distance, the same beam in the sonic crystal propagates without sensible broadening 
over much longer distances. Diffractive broadening in homogeneous medium is quantified by 
the Rayleigh distance Ld = ka'^/2, where a is the radius of the source, and corresponds to 
the distance from the source after which the beam broadens in a factor of V2. For the two 
nondiffractive frequencies evaluated above, the Rayleigh distances are 7.3 cm for the first 
case, and 10.3 cm for the second case. 

IV. ANALYTICS FOR NONDIFFRACTING BEAMS 

The precise analysis of an arbitrary field structure inside the crystal can only be performed 
by considering the field expansion into an infinite set of modes of the crystal (as stated by 
the Bloch theorem). The form given by Eqs. (0))-® must be assumed, whose unknown 
amplitudes can be numerically evaluated. This is the basics of the PWE method used in 
Sec. II for evaluate the band structure and dispersion curves of the crystal. However, it is 
possible to develop an analytical theory if we consider the particular case of a nondiffractive 
beam, since this nondiffractive solution appears mainly due to the coupling of three modes, 
the homogeneous one and the next low order modes. This situation is illustrated in Fig.5, 
where the three intersecting circles, corresponding to the spatial dispersion curves of the 
three modes (those with wavenumbers k, k + Gi and k + G2) give rise to the nondiffractive 
effect. Due to the interaction between the different spatial modes the degenerancy at the 
intersections of the spatial dispersion curves is lifted, and the fiat areas in the dispersion 



8 



curve can develop. The radiation belonging to these interacting modes is the most relevant 
for the deformation of the dispersion curves and to the appearance of the fiat segments, i.e. 
is responsible for the nondiffr active propagation. Therefore the other modes are irrelevant 
in the "nondiffr active" area (shadowed region in Fig. 5), and the Bloch expansion can be 
simplified as 

p(r)= (j9o+pie*^-^+P2e*^^-^), (10) 

where Gi and G2 are the basis vectors of the reciprocal space. 

Note that since the nondiffr active beam is expected to be highly directive, only G vectors 
directed to the same direction as the wavevector k are relevant in the analysis. The material 
parameters (being real functions) must be however expanded into at least five modes. In 
particular, the inverses of density and bulk modulus will be assumed of the form 

p{ry' = {bo + + 6+26*''^" + 6_ie-*^^" + fe.ae"*^^") , (11a) 

B{r)-^ = {bo + fo+ie^^^-"" + 6+26*^2-' + b^ie^'^' " + b.2e~'^^ ') , (lib) 

where the notation p±j = p±g^, with j=l,2 has been used. Substituting Eqs. (fTT|) in Eq. 
(ini), and collecting the terms at the same exponents (those with wavevectors k, k + Gi and 
k + G2), we obtain the following coupled equation system, 

= LJ^ {pobo + Pib-i + P2&-2) - k%po - k (k + Gi) pip_i - k (k + G2) P2P-2, (12a) 
= cu^ {p^bo + pob+i) - (k + Gi)' pipo - k (k + Gi) poP+i, (12b) 
= (P2&0 + Pob+2) - (k + G2)^ P2P0 - k (k + G2) PoP+2- (12c) 

Equations (jl2j) are still too complex to lead to analytical results. However, for small 
values of the filling fraction / the solid inclusions can be considered as a perturbation of the 
homogeneous fluid medium, and an assymptotic analysis near the bandgap is justified. Next 
we show that, in this limit, it is possible to obtain a simple relation between the frequency 
and wavenumber of the nondiffractive beam and the filling fraction / characterizing the 
crystal. 

First, note that in the case of small / (i.e. when tq << a) and materials with high- 
contrast, where ph « Pc and « (as occurs, e.g., for steel cylinders in water), the 
coefficients of the medium expansions in Eqs. ((Tj) and (jHl) simplify to po = &o = 1 — / and 
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Pi = hi = - 



/, for i = ±1, ±2. Then, Eqs. fll2p . expressed in matrix form, read 



(l-/)(^2-k2) /(k(k+Gi)-a.2) 
/ (k (k + Gi) - u') (1 - /) - (k + Gi)') 
^/(k(k+G2)-u;2) 



(l-/)(c.2-(k + G2f) J \P2 J 




0. 



(13) 



The aim is to obtain the values of u and k for which the beam propagates without diffrac- 
tion. For a crystal with square symmetry, the direction of the nondiffractive propagation 
with respect the crystal axes can be obtained from the isofrequency curves evaluated in Sec. 
II. For the first band (Figs. 3(a) and 4) nondiffractive propagation occurs for beams propa- 
gating at 45° with respect to the crystal axes, i.e. in the [1,1] direction. For our analysis, is 
convenient to consider the beam axis to be coincident with the z direction, and define a set 
of unitary basis vectors as Gi = (—1, 1)/V2, G2 = (1, l)/-\/2. In this way, all magnitudes 
in reciprocal space are normalized by vr/a. 

For small /, one also expects that the parameters corresponding to the nondiffractive 
regime take values close to those in the bandgap (near the corner of the first Brillouin zone; 
see Fig. 2). The wavenumber corresponding to the first bandgap is then Kb = (0, 1)7^2 
(remind that, in normalized reciprocal space, K = ka/27r). In order to investigate the 
behaviour of dispersion curves close to this point, we consider waves with wavevector K = 
Kb(1 — 6K.) with 6K. = {6Kx,SKz) representing small deviations. We further assume that 
the frequency is close to- (but less than) that corresponding to the bandgap, Q = Qb{^ — S^), 
with the normalized gap frequency given by fi^ = I/V2. 

The solvability condition of Eq. ()13|1 results from equating to zero the determinant of the 
matrix, and leads to the relation in the form F {6Q, 6Kz, SK^; f) = 0. Expanding for small 
5K = {6Kz, SK^), an analytical transversal dispersion relation 6Kz {SK^) is obtained, which 
allows to calculate the diffraction coefficient as the curvature of the transverse dispersion 
curve, i.e., D = {l/2)d'^5Kz/ dSK^. The nondiffractive point corresponds to D = 0. This 
expression is analytical but still cumbersome. However, assuming that / ~ (^{e^) and 
5Vt ~ 0{£), where e is a smallness parameter, the following simple analytical relation is 
obtained at the leading order: 



= + o if") . 



(14) 
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Also the wavenumber of the nondiffractive beam can be analytically evaluated. Substi- 
tuting Eq. ()14|) in the solvability condition of Eq. |T3j) . and assuming the above smallness 
conditions, it follows that, 

SK^ji^n = f^'-f^' + lf + Oif/'). (15) 

For the second band, a similar analysis can be performed. The three-mode expansion is 
illustrated in Fig. 6. In this case it is more convenient to use the vector basis Gi = (1, 0) and 
G2 = (0, 1), and now the nondiffractive beam propagation occurs in a direction coincident 
with one of the lattice vectors. The parameters of the gap are in this case = (1, 0) and 
= 1- An assymptotic analysis similar to that performed above for the first band, shows 
that = ^^TVD and 5K^]-, = 5Kj^^. Then, from the analytics follows that, under the 

limit of the weak modulation that / ~ Ole"^) and 6Q ~ C>{e), the zero diffraction point 
for both bands are situated equally, however with respect to the corresponding bandgap: 
the diffraction in the first band disappears just below the first bangap (by S^Ind = f^^^), 
and in the second band just below the second bandgap by the same value (by SQ^o)- The 
wavevector shifts are also equal for both cases. As a consequence, Eqs. fll4|) and ()15p are 
valid for both bands. 

These analytical predictions have been checked numerically. In Fig. 7 the analytical re- 
sults given by Eqs. (fT^ and are compared with the corresponding numerical results 
(with symbols) obtained with the PWE method using 361 modes. The curve labeled (a) 
corresponds to the normalized frequency shift, measured with respect to the bandgap, and 
the curve labeled (b) to the wavenumber shift, for zero diffraction point in the first band 
(circles) and the second band (squares). The simple expressions obtained under the three- 
mode expansions are in very with good agreement with the numerical results, even for the 
moderate (not very small) values of the filling factor /. 

V. CONCLUSIONS 

Concluding, we have demonstrated theoretically the possibility of nondiffractive propa- 
gation of acoustic beams through sonic crystals. We show the nondiffractive propagation 
for both propagation bands: for the first band, where the nondiffractive propagation occurs 
along the diagonals of the lattice, and for the second band, where diffraction disappears 
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along the principal vectors of the lattice. The diffraction disappears for frequencies just 
below the corresponding bandgaps, with equal frequency shift for both cases given by a 
universal and very simple expression: SflND = P^^. 

The other universal relation (jl5|) . for the shift of the wavenumber, which in simplified form 
reads Sk^D = f^^^, allows to evaluate the width of the nondiffractively propagating beams. 
Indeed the half-width of the platoe of spatial dispersion curve is roughly equal to (slightly 
less than) 5kND- This means that beams with normalized size d ^ 2T{/5kN£, ^ 27r/~^/^ 
can propagate nondiffractively over large distances (comparing with the diffraction length 
of the corresponding beam in the homogeneous materials). In the terms of non- normalized 
coordinates, the minimum size of the beam corresponds to d \/2a f^"^/^ . For a filling 
factor / = 0.114, corresponding to the numerical simulation in Fig. 4, the width of the 
nondiffractive beam predicted by this expression results d ^ 6a, i.e. to 6 spatial periods of 
the lattice. This result is in good agreement with the width of the beam observed in Fig. 4. 
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FIGURE CAPTIONS 



Fig.l. Geometrical interpretation of diffraction of waves propagating along the z axis: a) 
positive, or normal diffraction in propagation through homogeneous materials; b) negative, 
or anomalous diffraction; c) zero diffraction. The area of negligible diffraction (for evaluation 
of the minimum size of the nondiffractive beam) is indicated. 

Fig. 2. Band structure for steel cylinders in water, for r = 1 mm, a = 5.25 mm, as calcu- 
lated by the expansion into Bloch modes (4)- (7). The solid squares mark the nondiffractive 
points (see Fig. 3). 

Fig.3 Isofrequency lines, evaluated for a — 5.25 mm and r — 1 mm, for the first (a) 
and second (b) bands, centered at F point, as calculated by the expansion into Bloch modes 
(4)-(7). 

Fig. 4. Numerical FDTD simulation of the nondiffractive beam, for the first two bands. 
Upper row corresponds to the field radiated by the source without crystal, lower row to the 
nondiffractive propagation in the [1,1] (left) and [1,0] (right) directions, at the frequencies 
equal to f — 154 KHz and / = 217 KHz as predicted by the theory. Gray levels are in 
decibel scale, and the coordinates in meters. The other parameters are as in Figs. 2 and 3, 
i.e. steel cylinders in water are simulated with r — 1 mm and a — 5.25 mm. 

Fig. 5. Schematic picture showing the nondiffractive region (shadowed area) resulting 
from the interaction of three modes. The square represents the limits of the first Brillouin 
zone. The insert illustrates the lift of the degenerancy at the cross sections of the dispersion 
curves and the formation of the Bloch modes. The upper Bloch mode can develop flat 
segments depending on the interaction strength, as the degree of the lift of degenerancy is 
proportional to the interaction strength. 

Fig. 6. Schematic picture showing the nondiffractive region (shadowed area) in the second 
propagation band. Everything as in Fig. 5. Here the most relevant modes are k + Gi and 
k±G2. 

Fig. 7. Dependence of the frequency (a) and wavenumber of nondiffractive beam, mea- 
sured with respect to the bandgap values, as results from numerical calculation (symbols) 
and the analytical expressions given in Eqs. (12) and (13). The open circles and squares 
correspond to the parameter values used for FDTD calculation of nondiffractive propagation 
(Fig. 4) in the first and second band respectively. 
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